//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "statics.hpp"
#include "screen.hpp"
#include "dynamics.hpp"
#include "consts.hpp"

// test for e field along axis of a circular charge

#define RADIUS 1.0
#define CHARGE 0.00000001
#define DIST0 1e-5
#define DIST1 0.1
#define DIST2 0.65
#define DIST3 3.65
#define DIST4 8.65
#define ALLOWED_FACTOR 1000.0 // allow difference by 1/1000 worth

static void test_against (const char *where, const vect3d& ef, const vect3d& t_ef)
{
  prec_t diff = (ef - t_ef).length();
  if ( diff > t_ef.length() / ALLOWED_FACTOR )
    {
      printf ("ERROR_%s=fields differ too much\n", where);
    }
  else
    {
      printf ("OK_%s=fields match up\n", where);
    }
  char buf[128];
  printf ("EF_%s=%s\n", where, ef.sprint (buf));
  printf ("TEF_%s=%s\n", where, t_ef.sprint (buf));
  printf ("MARGIN_%s=%g (%g%%)\n", where, prec2double(diff), prec2double(100.0 * diff / t_ef.length()));
}

// from http://hyperphysics.phy-astr.gsu.edu/hbase/electric/elelin.html#c2
// Ez = kQz/r^3

#define KE (MAGNETIC_CONSTANT*SPEED_OF_LIGHT*SPEED_OF_LIGHT)

static prec_t Ez (prec_t d)
{
  prec_t r2 = d*d + RADIUS*RADIUS;
  prec_t r3 = r2 * prec_t_sqrt(r2);
  return KE * CHARGE * d / r3;
}

int main (int argc, char *argv[])
{
  Statics statics;
  ElectroDynamics dyn(statics);
  StaticRing *pc = new StaticRing(0, vect3d(0, 0, 0), RADIUS, 0, vect3d(0.0, 0.0, 1));
  pc->setCharge (CHARGE / (PREC_PI * 2 * RADIUS));
  statics.addStaticElement (pc);

  vect3d ef;
  dyn.getEField (vect3d(0, 0, DIST0), ef);
  test_against ("DIST0", ef, vect3d (0, 0, Ez(DIST0)));
  dyn.getEField (vect3d(0, 0, DIST1), ef);
  test_against ("DIST1", ef, vect3d (0, 0, Ez(DIST1)));
  dyn.getEField (vect3d(0, 0, DIST2), ef);
  test_against ("DIST2", ef, vect3d (0, 0, Ez(DIST2)));
  dyn.getEField (vect3d(0, 0, DIST3), ef);
  test_against ("DIST3", ef, vect3d (0, 0, Ez(DIST3)));
  dyn.getEField (vect3d(0, 0, DIST4), ef);
  test_against ("DIST4", ef, vect3d (0, 0, Ez(DIST4)));

  return 0;
}

